Trap spaces of multi-valued networks: definition, computation, and applications

Abstract Motivation Boolean networks are simple but efficient mathematical formalism for modelling complex biological systems. However, having only two levels of activation is sometimes not enough to fully capture the dynamics of real-world biological systems. Hence, the need for multi-valued networks (MVNs), a generalization of Boolean networks. Despite the importance of MVNs for modelling biological systems, only limited progress has been made on developing theories, analysis methods, and tools that can support them. In particular, the recent use of trap spaces in Boolean networks made a great impact on the field of systems biology, but there has been no similar concept defined and studied for MVNs to date. Results In this work, we generalize the concept of trap spaces in Boolean networks to that in MVNs. We then develop the theory and the analysis methods for trap spaces in MVNs. In particular, we implement all proposed methods in a Python package called trapmvn. Not only showing the applicability of our approach via a realistic case study, we also evaluate the time efficiency of the method on a large collection of real-world models. The experimental results confirm the time efficiency, which we believe enables more accurate analysis on larger and more complex multi-valued models. Availability and implementation Source code and data are freely available at https://github.com/giang-trinh/trap-mvn.


Introduction
Boolean networks are simple but efficient mathematical formalism for modelling, analysing, and controlling complex biological systems (Schwab et al. 2020). Beyond systems biology, they have been widely applied in various areas from science to engineering (Schwab et al. 2020). Boolean network models of biological systems represent genes (or other species) as nodes that can take Boolean values: 1 (active) and 0 (inactive). However, having only two levels of activation may not be enough to fully capture the dynamics of real-world biological systems (Schaub et al. 2007). There are many examples (Schaub et al. 2007;Didier et al. 2011;Mushthofa et al. 2018) where the dynamics of the system can only be modelled by considering more than two activation levels. Hence, there is a crucial need to study multi-valued networks (MVNs), which are a generalization of Boolean networks (Naldi et al. 2007;Schaub et al. 2007).

Related work
Despite the importance of MVNs, only limited progress has been made on developing theories, analysis methods, and tools that can support them (Mushthofa et al. 2018). First, besides simulation, the analysis of logical models is mostly based around 'attractor' computation, since those correspond roughly to observable biological phenotypes (Schwab et al. 2020). For example, in gene regulatory and signalling networks, attractors can correspond to cell types, cell fates, and cyclic behaviour (e.g. circadian rhythms and cell cycles). Hence, analysis of attractors could provide new insights into systems biology. However, finding all attractors of a logical model (even for the Boolean case) is challenging due to the complex dynamics of models (Schwab et al. 2020). The recent study of trap spaces of Boolean networks (Klarner et al. 2017) made a real breakthrough in the field of systems biology, as minimal trap spaces provide very good approximations of attractors and are much easier to compute. However, there has been no similar concept defined and studied for MVNs to date.
Furthermore, other biological properties, such as the gradual commitment of a cell to a specific phenotype, can be revealed through the lens of 'succession diagrams' constructed from the networks' 'maximal' trap spaces (Rozum et al. 2021(Rozum et al. , 2022. As such, minimal trap spaces are not the only relevant form of a trap space. In the Supplementary Data (Section S4.2), we give a more detailed discussion of trap space applications in biological modelling.
One notable issue is that the current supporting methods for MVNs cannot handle large and complex models (Naldi et al. 2007 ;Mushthofa et al. 2018). This issue also prevents the modellers from building such models (Naldi et al. 2007), which could provide more accurate insights. Finally, a popular research direction is to convert an MVN to a Boolean network with similar dynamical behaviour, then applying the rich set of analysis methods/tools designed for Boolean networks. However, the existing Boolean encoding methods (e.g. the Van Ham Boolean mapping and its variants;Ham 1979;Didier et al. 2011) may not cover the full set of dynamics of the original MVN (Didier et al. 2011), and the encoding may even hinder the efficiency of the Boolean network methods/ tools (Mushthofa et al. 2018). It is also worth noting that all mentioned encoding methods support only unitary MVNs. We believe that it is possible to develop direct and efficient methods for MVNs.

Our contributions
In this work, we study both general and unitary MVNs. Note that this inclusion of general networks is important not only in theory, but also in many biological applications (see Supplementary Section S1). First, we generalize the concept of trap spaces in Boolean networks to that in MVNs. Second, we prove several properties of trap spaces in MVNs including (i) the separation of minimal trap spaces, (ii) their relation with respect to the Van Ham Boolean mapping, and (iii) the characterization of trap spaces through Petri net (PN) siphons. Based on this characterization, we propose a new method utilizing answer set programming (ASP; Gebser et al. 2011) for computing different types of trap spaces of MVNs, including generic trap spaces, maximal trap spaces, minimal trap spaces, and fixed-points (a special sub-type of trap spaces). For fixed-points, we also consider another new method, which relies on the characterization of deadlocks of the PN encoding (Liu and Barkaoui 2016). We evaluate the method on a collection of real-world MVNs and show its applicability in treatment prediction on a case study of Myc-deregulation in breast cancer.

Paper outline
This article is structured as follows: In Section 2, we formalize the notion of trap spaces in MVNs and their relationship to the PN siphons. We propose a method based on answer-set programming for computing such trap spaces efficiently. In Section 3, we present a Python package trapmvn that implements this method for both SBML-qual and BMA models. We demonstrate the practical performance of the method on a wide selection of real-world models. We then use trapmvn as part of a case study to identify viable intervention targets in a model of Myc-associated deregulation in breast cancer. Finally, Section 4 discusses the scalability of the method, its future prospects, and highlights the role of trap spaces in reliable long-term behavioural analysis of logical models.
This article also has an associated supplement, in which we include the non-essential technical details of our methodology, as well as the full results of our benchmarks and case study. This article refers to this supplement where appropriate.

Methods
In this section, we define trap spaces of MVNs, discuss some of their theoretical properties, and finally present a method based on ASP for computing trap spaces of symbolically represented MVNs.

MVNs
DEFINITION 1. An MVN is a triple M ¼ ðV; K; FÞ such that: . . . ; v n g is an indexed set of nodes (variables). • K ¼ fK 1 ; . . . ; K n g is an indexed set of integer intervals, representing the domains of variables v i 2 V. • F ¼ ff 1 ; . . . ; f n g is an indexed set of update functions of variables v i 2 V. Each f i has a signature f i : P n j¼1 K j ! K i .
With a slight abuse of notation, we write v i to denote both the network node and the associated integer variable. Furthermore, we can also write K vi and f vi to denote the domain and update function of v i , respectively. We define x 2 Q n j¼1 K j as the 'state' of the MVN, with x i (or x vi ) denoting the value of v i in the state x. We can also write x as a vector ½x 1 ; . . . ; x n . We then write S M ¼ Q n j¼1 K j to denote the set of all states ('state space') of network M. Finally, note that there are different possible formulations of update functions F, including fuzzy logic (Mushthofa et al. 2018), rule-based descriptions (Naldi et al. 2007; Delaplace and Ivanov 2020), or arithmetic expressions (Benque et al. 2012). We will return to this aspect when we discuss the encoding and manipulation of networks by our method. The dynamics of an MVN are studied through its 'statetransition graph', STG M ¼ ðS M ; !Þ. In particular, a 'trap' set of STG M is a set of states that is closed with respect to !. A trap set is called an 'attractor' if it is minimal, i.e. there is no other trap set that is a proper subset of this set. However, there are different ways of defining STG M , possibly leading to significantly different behavioural features.
First, we divide MVNs based on the variable update scheme as 'general' (Mushthofa et al. 2018) and 'unitary' (Schaub et al. 2007;Delaplace and Ivanov 2020). In the general networks, the variable changes follow the update functions exactly. Meanwhile, in the unitary networks, the value of v i can only change by 'one level' at a time (i.e. þ1 or -1). It has been argued that unitary networks better capture the continuity of biological interactions (Schaub et al. 2007), but is not true universally (Sun et al. 2014).
Note that the stepwise notion of unitary update can be in fact encoded directly into the update functions F. Specifically, we definef v to represent the update function which adheres to the chosen variable update. That is, The second aspect of MVN semantics is the concurrency of updates. Here, two prevalent approaches are the 'synchronous' (all variables update together) and 'asynchronous' update (exactly one variable updates non-deterministically in each step).
Note that for the same MVN, different update schemes may result in vastly different state-transition graphs (Thomas 1991). Consider the MVN M with V ¼ fv 1 ; v 2 g, K 1 ¼ f0; 1g, and K 2 ¼ f0; 1; 2g (Example 2 in the Supplementary Data). Figure 1

Trap spaces of MVNs
Overall, this variability under different update schemes motivates the study of 'trap spaces'. In the Boolean case, trap spaces have been shown to be a good approximation of network attractors, regardless of the chosen update concurrency (Paulevé et al. 2020). Here, we thus establish trap spaces for MVNs with a similar goal in mind: DEFINITION 2. A space m of S M is a mapping which assigns each v i 2 V a non-empty subset of K i : mðv i Þ K i . We write S ? M to denote the set of all spaces of M.
With a slight abuse of notation, we can also interpret m as a subset of S M , writing x 2 m for a state x 2 S M when x v 2 mðvÞ for all v 2 V. Subsequently, we can define trap spaces as follows: DEFINITION 3. A trap space is a space m such that for every state x 2 m and variable v 2 V, we havê f v ðx v Þ 2 mðvÞ.
In other words, in all states represented by m, all update functions only produce values that stay within m, meaning no STG M (regardless of the update concurrency) can contain a transition which leaves m. Consequently, trap spaces of an MVN are independent on its concurrency update, is not true for the case of attractors. Note the use off v , meaning that the notion of trap space respects the choice of unitary or general variable update.
When interpreting spaces m 2 S ? M as subsets of S M , we can trivially establish a sub-space relation on S ? M , which allows us to reason about 'minimal' and 'maximal' trap spaces: Note that for the notion of maximal subspace, we require that m & S M . Otherwise, a trivial space m ? which represents all network states (i.e. m ? ðvÞ ¼ K v ) is always a trap space, and always a superset of all other trap spaces.
For illustration, reconsider the example MVN ( Fig. 1). Its general counterpart has six trap spaces as follows: Out of these, three are minimal trap spaces (m 1 , m 2 , and m 3 ) and two are maximal trap spaces (m 4 and m 5 ). Meanwhile, its unitary counterpart has two minimal trap spaces (f0gf1g and f1gf1g) and three maximal trap spaces (f0; 1gf0; 1g, f0; 1gf1; 2g, and f0gf0; 1; 2g). Finally, let us observe the following: THEOREM 1 (Separation of minimal trap spaces). Let m 1 and m 2 be two distinct minimal trap spaces of M. Then m 1 \ m 2 ¼ 1.
The proof is available in the Supplementary Data as Theorem 3.
By Definition 3 and the definition of a trap set, a trap space is also a trap set of an MVN for any chosen update concurrency. As a consequence, each trap space must contain at least one attractor of the MVN. By Theorem 1, any two distinct minimal trap spaces are disjoint; thus, the attractors contained in them are also disjoint. Hence, we can conclude that regardless of the chosen update concurrency, minimal trap spaces can be used as approximations of attractors in MVNs: The number of minimal trap spaces under-approximates the number of attractors (there can be other attractors that are not contained in minimal trap spaces), and each minimal trap space over-approximates at least a single attractor. Supplementary Section S3.3 presents additional examples of relations between attractors and minimal trap spaces.

Translation to Boolean networks
When K i ¼ f0; 1g for all v i 2 V, we can refer to M as a 'Boolean network' N ¼ ðV; FÞ instead. For Boolean networks, the notion of trap spaces is well known, including efficient tools for computation of minimal and maximal trap spaces (Klarner et al. 2017;Paulevé et al. 2020;Trinh et al. 2022).
There are several works (Didier et al. 2011;Delaplace and Ivanov 2020) attempting to encode the multi-valued dynamics of an MVN into a Boolean network. Probably, the most widely used being the Van Ham encoding (Ham 1979), implemented in the tools GINsim (Chaouiya et al. 2011a(Chaouiya et al. , 2011b and bioLQM (Chaouiya et al. 2013). The core idea of the Van Ham encoding is to expand each multi-valued variable v 2 V into jK v j À 1 Boolean variables, such that an integer value v ¼ k is encoded as the truth value of the first k Boolean variables. The main advantage of this encoding is that any unitary state change only involves a single Boolean variable. The disadvantage is that the encoding admits invalid states which do not correctly encode any integer and may interfere with the actual dynamics of the model. More details about the Van Ham encoding are given in the Supplementary Section S2.3.
It is known that the Van Ham encoding preserves attractors of unitary networks (Didier et al. 2011). Nevertheless, it is unclear whether this mapping also preserves trap spaces in general. To answer this question, we present the following This proposition states that the Van Ham encoding is still applicable for computation of minimal trap spaces of unitary networks, but cannot be used beyond this particular problem class.

Computing trap spaces through PN siphons
To tackle the problem of trap space computation of MVNs, we instead propose a different approach, based on an encoding of the problem into an answer-set programming query.

PN encoding of MVNs
First, let us establish a characterization of trap spaces through conflict-free PN siphons.
DEFINITION 5. A one-safe PN is a bipartite-directed graph P ¼ ðP; T; WÞ, where P and T are disjoint finite sets of vertices called places and transitions, respectively. Set W describes the arcs between places and transitions: W ðP Â TÞ [ ðT Â PÞ.
A marking M of a one-safe PN is a subset M P. For x 2 P [ T, we write pred(x) and succ(x) to denote the predecessors and successors of x with respect to W (this notation naturally extends to subsets of P [ T). The dynamics of a PN are dictated by the 'firing of transitions', such that a transition t 2 T can be fired non-deterministically in the marking M if predðtÞ M. The result is a new marking M 0 ¼ ðM n predðtÞÞ [ succðtÞ. This process defines a statetransition graph with vertices as the possible markings.
In the appendix of Chatain et al. (2014), the authors present a PN encoding of the 'asynchronous' MVN dynamics, which we briefly recall here. Note that this is not the only PN encoding of MVNs (see, e.g. Chaouiya et al. 2004Chaouiya et al. , 2011aChaouiya et al. , 2011b, but to the best of our knowledge, the other encodings are not particularly suitable for the characterization of trap spaces. Let P M denote the one-safe PN encoding of an MVN M based on Chatain et al. (2014). The places P of P M contain one place for every level of every MVN variable: P ¼ [ v2V fp v¼i ji 2 K v g. The set of transitions and arcs is then constructed such that the firing of a transition mirrors a possible update of an MVN variable. For example, let u; v 2 V be such that for every x 2 S M , we have We then create a transition t with predðtÞ ¼ fp v¼k ; p u¼l g and succðtÞ ¼ fp v¼o ; p u¼l g. Such transition 'moves' a token from p v¼k to p v¼o under the assumption u ¼ l. Note that there are many possible combinations of sets T and W that are valid for a particular M. Later, we show how we compute these sets in our case.
Finally, let m 2 S ? M be a space of M. We write that a set of places M½m P is the 'mirror' of m when for every v 2 V and i 2 K v ; i 2 mðvÞ () p v¼i 6 2 S. That is, M½m contains the places corresponding to the 'inverse' of m. Observe that a state x 2 S M is also a (trivial) space, and hence has a mirror M½x.

Siphon characterization of MVN trap spaces
DEFINITION 6. A PN siphon S P is a set of places such that for all t 2 T; ðS \ succðtÞÞ 6 ¼ 1 implies ðS \ predðtÞÞ 6 ¼ 1.
Intuitively, a siphon S preserves the condition that if all places in S are unmarked (M \ S ¼ 1), the siphon remains unmarked. Furthermore, for a PN P M , we say that a siphon is 'conflict-free' when for all v 2 V, we have that fp v¼i ji 2 K v g n S 6 ¼ 1. That is, a conflict-free siphon must not contain 'all' places encoding a particular variable. Intuitively, every conflict-free siphon of P M represents a mirror of 'some' space m 2 S ? M (however, not all mirrors of spaces are siphons).
With this knowledge, we can observe the following theorem: THEOREM 2. Let M be an MVN and P M its PN encoding.
Then a space m 2 S ? M is a trap space if and only if its mirror M½m is a conflict-free siphon of P M .
The proof is available in the Supplementary Data as Theorem 5.
Finally, note that we can define a partial order on siphons based on the subset relation, just as we did for trap spaces. Based on the definition of a mirror and the correspondence we just proved, we can easily deduce the following (see the proofs of Theorems 6 and 7 in the Supplementary Data):

PROPOSITION 2. A trap space m of an MVN M is minimal if
M½m is a maximal conflict-free siphon of P M . A trap space m is maximal if M½m is a minimal conflict-free siphon.
Observe that there are existing methods based on ASP (Trinh et al. 2022) and SAT solving (Nabli et al., 2016) that can handle enumeration of maximal (respectively, minimal) siphons and have been successfully used on Boolean networks before (Trinh et al. 2022). Our result shows that these methods should be also applicable to MVNs. Specifically, we aim to extend the ASP encoding from (Trinh et al. 2022) in order to support both maximal and minimal trap spaces of both general and unitary MVNs.

Siphon computation through ASP encoding
For now, we neglected the question of update function representation for MVNs, as this depends strongly on the chosen input format. In our work, we consider both SBML-qual (Chaouiya et al. 2013) and BMA (BioModelAnalyzer) (Benque et al. 2012) formats. These differ significantly in terms of both format and capabilities.
Ultimately, we represent an update function f i using a series of 'binary decision diagrams' (BDDs) (Bryant 1986) B 1 ; . . . ; B k . The network variables are encoded into Boolean variables [ v2V fp v¼i ji 2 K v g, just as in the case of the PN encoding. Subsequently, each BDD in the series gives the necessary and sufficient conditions for achieving a particular result level: Here, K i ¼ ½y 1 ; . . . ; y k and B i are individual BDDs, such that each BDD can be interpreted as a Boolean function on the encoded state x. We also require that for each state, there is exactly one BDD B j s.t. B j ðxÞ ¼ 1, hence the whole function is well defined. Also, note that such f i can be easily transformed intof i , respecting either general or unitary variable update as desired.
This representation is relatively straightforward to obtain for SBML-qual models, as each update function is given as a list of Boolean terms over standard equality/inequality integer propositions. In this representation, the output values can repeat across multiple terms, and the function is not required to be exhaustive (there is a 'default' output value assigned to remaining inputs). However, this is easy to amend once each term is transformed into a symbolic BDD representation (see Supplementary Data).
For BMA, the translation process is more involved, as update functions are described through a language of algebraic expressions including addition (þ), subtraction (-), multiplication (Á), and division (=), as well as other special functions like average and rounding. Furthermore, BMA employs a normalizing transformation on function inputs when the input domain differs from the output domain. That is, an input variable v in the range K v ¼ ½x; y is normalized to the range K u ¼ ½a; b when used in the update function f u . As such, while BMA only admits integer variables, the update functions are more akin to rational functions.
While there are frameworks which partially support symbolic evaluation of such functions, e.g. algebraic decision diagrams (Bahar et al. 1997), we are not aware of any implementation that would support all the operations required by BMA. We thus opted to enumerate the whole function table and re-encode it back into individual BDDs.
Finally, to encode the dynamics off into a PN, we consider all variable updates ðv ¼ lÞ ! ðv ¼ kÞ (with k being l61 when M is unitary), and then enumerate all satisfying partial valuations w of the BDD p v¼l^Bk (B k being the BDD off for output level k). For every valuation w, a transition t 2 T is created which moves the token from p v¼l to p v¼k while ensuring that for every other relevant u 2 V, place p u¼wðuÞ contains a token.
Note that the number of satisfying partial valuations w of a BDD (and transitively, the number of PN transitions) depends on the 'ordering' of Boolean variables within the BDD. However, computing the optimal variable ordering is a known non-trivial problem. To reduce the number of transitions in the resulting PN, we always test k randomized orderings for every update function and pick the one which produces the most compact PN.
Once the PN encoding P M is completed, the method proceeds based on the ASP encoding of the PN siphon problem proposed by Trinh et al. (2022). This encoding produces a query which can be processed by an ASP solver such as clingo (Gebser et al. 2011), enumerating all maximal/minimal siphons of P M . Technical aspects of this process are given in the Supplementary Data.

Results
We now present the results of our computational experiments. First, we establish the benchmark model dataset and describe the overall implementation of our tool trapmvn. We then present the performance evaluation of trapmvn in relation to the benchmark models and other tools (where applicable). Finally, we use the trap spaces arising under a wide range of therapeutic interventions in a Myc-deregulation model of breast cancer to assess the viability and reliability of such interventions.

Experiment setup
In our testing, we utilize SBML models from the BBM benchmark (https://github.com/sybila/biodivine-boolean-models/; note that BBM primarily publishes Booleanized SBML models, but the original multi-valued SBML files are also available.), which includes models from the GINsim tool (Chaouiya et al. 2011a(Chaouiya et al. , 2011b repository as well as other independently sourced models. Furthermore, we include models from the BMA tool (http://biomodelanalyzer.org) repository (Benque et al. 2012). Disregarding trivial cases, we are left with 26 benchmark models. To the best of our knowledge, this is a highly representative sample of multi-valued models currently available in literature. The details of each benchmark model, technical description of the performance testing, as well as full results are given in Section 6 of the Supplementary Data.

The trapmvn package
We implement our method as a stand-alone open-source Python package trapmvn. The package provides a basic parser for both SBML and BMA models (JSON or XML). It implements the symbolic encoding for both formats using the BDD data structure from the AEON.py package (Bene s et al. 2022). From this symbolic encoding, we build either a 'general' or a 'unitary' one-safe PN. Such PN can be then encoded into the ASP query for the respective problem class and processed by the solver clingo (Gebser et al. 2011). All these steps are available to the user through several Python classes, but can be also handled directly by a simple command line interface. Finally, we also support export of the symbolic model representation back into SBML. This allows us to convert BMA models to SBML (to the best of our knowledge, this is not supported by any other tool).

Performance evaluation
Due to the lack of both theory and available tools, we are not aware of any existing study of trap spaces in MVNs. Nevertheless, as we have shown, the Van Ham encoding (Ham 1979) preserves the 'minimal' trap spaces under unitary semantics. Tools for computing minimal trap spaces of plain Boolean networks can be thus used in this case. We use bioLQM (Chaouiya et al. 2013) to compute the Van Ham encoding when necessary.
We can also consider fixed-points as a special (simpler) class of trap spaces. We can thus also compare with tools that specialize in this type of problem. Consequently, our Trap spaces of multi-valued networks i517 performance evaluation consists of three parts: 'minimal' and 'maximal' trap spaces, and 'fixed-points'.

Minimal trap spaces
To test the performance of trapmvn when computing minimal trap spaces, we compare its runtime to two state-of-theart tools for trap space computation in Boolean networks, trappist (Trinh et al. 2022) and mpbn (Paulevé et al. 2020) (using the Van Ham encoding). Because the number of trap spaces can be large and the knowledge of all trap spaces is not always required, we consider two experiment settings: computing just 'one' and computing 'all' minimal trap spaces. The full results of this analysis are available in Supplementary  Tables S5 and S6. In general, trapmvn can handle all models in reasonable time for both the general and unitary semantics. Focusing on the unitary case, the difference in runtime when computing all minimal trap spaces is summarized in Fig. 2. Here, we see that trapmvn performs substantially better than mpbn, and is also always faster than trappist. However, the exact speed-up is hard to assess due to the presence of logarithmic time scales. As such, we also provide a simple box plot comparing the actual speed-up in Fig. 3. Here, we see that trapmvn outperforms both methods by a significant margin in both scenarios (i.e. first result and all results). We also compute the average speed-up weighted by the absolute runtime of each benchmark (i.e. longer running benchmarks are given higher weight), showing that when computing 'all' trap spaces, trapmvn is 3Â faster than trappist, and 58Â faster than mpbn.

Maximal trap spaces
The Van Ham encoding is not suitable for preserving 'maximal' trap spaces, and we are not aware of any existing method that can efficiently transform this problem to the domain of Boolean networks. As such, we cannot compare trapmvn to any other tool for this problem class. Nevertheless, we show that trapmvn can easily compute all maximal trap spaces of the considered benchmark models (Supplementary Table S7).

Fixed-points
Finally, for the case of fixed-points, there are two viable ASP encodings, one based on PN siphons and the other based on PN deadlocks. The encoding based on deadlocks is generally more efficient for this simpler problem class. As such, we implement both variants in trapmvn and compare the results to the Booleanized results from the tools trappist (Trinh et al. 2022) and mpbn (Paulevé et al. 2020), and multi-valued results from AN-ASP (Abdallah et al. 2017). Overall, our evaluation confirms the benefits of the deadlock-based encoding, showing a 2:6Â speed-up compared with the siphon encoding (Supplementary Table S8). Furthermore, we show that trapmvn also performs better than trappist or mpbn on fixed-point computation, but is outperformed by AN-ASP (1:3Â speed-up), since AN-ASP is optimized solely for this specific type of problem.

Therapeutic interventions of Myc-deregulation
To demonstrate the practical utility of trapmvn in biological modelling, we present a case study expanding on the findings of Kreuzaler et al. (2019). The authors of Kreuzaler et al. (2019) use BMA (Benque et al. 2012) to explore viable therapeutic interventions in a large computational model of breast cancer through the lens of network attractors. Here, we show how to more reliably interpret the model by focusing on trap spaces instead.

Modelling Myc heterogeneity in breast cancer
The Myc transcription factor is one of the key coordinators in cell proliferation and regeneration (Kortlever et al. 2017). As such, oncogenic deregulations of Myc are commonplace in many cancers, breast cancer in particular (Vita and Henriksson 2006).
Still, most tumours have been shown to consist of several genetically distinct mutants, only some of which exhibit Myc overexpression (Gerlinger et al. 2012;Heselmeyer-Haddad et al. 2012). Such heterogeneity can impede some treatments, but it can also enable new therapies that target the cooperation between the mutants (Marusyk et al. 2014;Kreuzaler et al. 2019).
In the case of Myc-related mutations, an overexpression of Myc is linked to super-competitive behaviour that causes the cancerous cells to outproliferate their healthy neighbours. However, the same overexpression is also linked to greatly increased predisposition to apoptosis. In Kreuzaler et al. (2019), the authors reveal a mechanism by which Myc high mutants survive using a supply of Wnt1 transcription factor produced by a different, Myc low mutant.
This process is demonstrated both experimentally on in vivo mouse models, as well as in silico on a large-scale   Trinh et al.
multi-valued computational model. Based on the intervention response observed in silico, the authors identify a viable therapy targeting COX2 and MEK transcription factors and validate this therapy in vivo.
The model itself consists of 72 variables, ranging from 4 to 7 levels (MYC IN-VIVO in Supplementary Table S5). It is based on known literature, with additional validation and tuning using several independent datasets. There are five model variants: A healthy 'wild-type' (WT) model, Myc low and Myc high variants where the tumour consists homogeneously of a single mutant, and finally, mix-Myc low and mix-Myc high variants, which describe a heterogeneous tumour. In this case, the interaction is given as an outside assumption; there is no single model consisting of both Myc low and Myc high mutants sharing a state space.
To study the effects of possible therapeutic interventions, the authors compare the 'synchronous' attractors approximated by BMA across a range of single and dual variable knockouts. The effect of each intervention on the real-world phenotypes is assessed through variables Apoptosis and Proliferation which are directly embedded into the network. While this methodology is viable, it has shortcomings which we hope to address.
First, while BMA should be capable of computing the 'exact' synchronous attractors, the results in Kreuzaler et al. (2019) are only based on an 'approximate' method, due to the (lack of) scalability of the exact method. Second, the 'synchronous' update scheme can miss plausible model behaviour due to artificial synchronization between variables (Schwab et al. 2020). Meanwhile, trap spaces are universal regardless of the chosen update scheme (Paulevé et al. 2020). Finally, for nontrivial attractors, the case study in Kreuzaler et al. (2019) only considers the 'average' values of Apoptosis and Proliferation, which can be a poor approximation of the model's actual admissible behaviour. This is despite the fact that BMA 'can' also compute the estimated attractor intervals.

Single intervention effects
We start by replicating the single-variable knockout interventions performed in Kreuzaler et al. (2019), but in the context of trap spaces. Highlights from this analysis are presented in Table 1 (full data are available in the Supplementary Data). These results are in many aspects comparable to the original results obtained through BMA, however, they often paint a more complete picture of the model's behaviour.
First, there are four instances where our 'exact' method actually improves the precision of the original 'approximate' result (marked with an asterisk). Note that the more specific result is in all cases at the edge of the interval established by BMA, which means that it substantially differs from the 'average' considered in Kreuzaler et al. (2019).
Second, there are clearly many cases where the intervention causes the appearance of a non-trivial trap space, signified by an interval instead of a fixed value. Knowledge of these intervals is crucial when interpreting the effectiveness of interventions.
For example, consider the value of Apoptosis for the knockouts of Mcl1 and COX2 in mix-Myc low . The average value is the same (i.e. 3), but the admissible interval is ½2; 4 for Mcl1 and ½0; 6 for COX2. The Mcl1 intervention guarantees Apoptosis ! 2, but we have no such expectation for COX2: Even though the 'best case' outcome is higher (i.e. 6 > 4), the system is not guaranteed to visit these high-value states sufficiently often to trigger apoptosis.

Dual interventions with reliable and opportunistic effects
To systematically rank the high number of possible dual interventions, we propose to score each intervention I with two metrics, 'reliability' and 'opportunity', denoted by rel(I) and opp(I), respectively. The intuition is that these scores should represent the 'worst' and 'best' case scenarios (in terms of Apoptosis and Proliferation) admitted by the model for each intervention, regardless of the considered update Trap spaces of multi-valued networks i519 concurrency. The exact definitions of these scores are given in Supplementary Section S5.4. We then focus on the admissible dual-knockout interventions: both knockouts must be druggable (according to Kreuzaler et al. 2019), and their combination cannot increase Apoptosis beyond level three in the WT model. This leaves 995 interventions for which we compute both scores. For presentation purposes, we sort the interventions by the 'average' of these two scores.
The eight best and the eight worst interventions are shown in Fig. 4. As we can see, some of the best interventions admit a non-trivial trap space for the mix-Myc low variant. However, even in this case, the reliability and opportunity scores are not vastly different. This raises a natural question regarding the prevalence of differences between reliability and opportunity scores in general. We further study this question in the Supplementary Data (Supplementary Section S5.4 and Supplementary Fig. S6, in particular), where we show that substantial differences between the two scores are in fact common even within the best scoring interventions.

Discussion
In this article, we formalized the notion of trap spaces in MVNs, then explored and proved properties of such trap spaces with applications in the analysis and control of MVNs. One notable property is that trap spaces of MVNs are independent to concurrency update schemes. We argued that, akin to Boolean networks, trap spaces can serve as approximation of network attractors. Furthermore, we showed that not all useful properties of trap spaces are preserved through a Boolean encoding of an MVN. For example, maximal trap spaces of unitary networks, crucial for the construction of the network's succession diagram (Rozum et al. 2021), are not preserved. Next, we made a connection between trap spaces of an MVN and siphons of its PN encoding. Based on this relationship, we proposed and implemented a new ASP method for computing different types of trap spaces of an MVN.

Method performance
We have evaluated the time efficiency of our method on realworld models collected from the literature. We show that the method scales well with the network size and it can handle large-scale realistic models for both the general and unitary semantics. The indirect approach (i.e. through a Boolean encoding) is only applicable for the case of fixed-points and minimal unitary trap spaces. In these cases, the direct approach (i.e. our method) outperforms the best indirect method.
In particular, we discuss the factors that contribute to the running time of each compared method (all are ASP-based) with respect to the minimal trap space computation. Aside from the absolute number of solutions, the practical complexity of an ASP query is affected by its number of atoms and its 'density', i.e. the ratio between the number of ASP rules and atoms. Through a systematic analysis (details are available in Supplementary Section S6.1), we see that an increase in density is always accompanied by increased runtime. Query density is clearly not the only indicator of ASP problem complexity. However, assuming we control for other relevant factors (solution count, tool/method, update scheme, model format, etc.), query density appears to be a relevant metric for comparing the complexity of two minimal trap space computation problems.
We also analyse the runtime of our method trap-mvn in the general and unitary cases. In general, trap-mvn needs more time for the general case than for the unitary case, as the PN encoding of a general MVN has more transitions than that of its unitary counterpart. The detailed discussions are given at the end of Section S6.1.
Finally, we have tried to also compare the performance of trap-mvn to BMA. However, in its current state, we were only able to run the 'approximate' attractor detection method, which is (as expected), much faster than any of the tested 'exact' methods. There appears to be no public documentation or reproducible artefact concerning the exact solver-based approach for the currently available version of BMA. However, the authors of Kreuzaler et al. (2019) note that the exact method was not able to efficiently analyse the Myc heterogeneity model. As such, it seems unlikely that it would provide competitive performance on the remaining models in our benchmark set.

Reliable identification of network interventions
Subsequently, we studied the practical applicability of our method on a model of Myc-deregulation from Kreuzaler et al. (2019). First, we indeed found instances where our exact method can improve the approximate results obtained in Kreuzaler et al. (2019). Furthermore, as we show in Table 1, the behaviour of the model can often admit significant fluctuations and uncertainty. These are not reflected in the average value within attractor states, as considered in Kreuzaler et al. (2019). This highlights the need for rigorous and exhaustive analysis of such behaviour uncertainty.
With this goal in mind, we propose a 'reliability' and 'opportunity' score to assess the viability of therapeutic interventions. Based on these scores, we see that even interventions that introduce some amount of uncertainty can be safely considered among the most viable (Fig. 4). Furthermore, we study how prevalent is the variability of intervention scores in our  Trinh et al.